clear
rng default

%% Load useful matrices
load G_round
x=1;
all_years=1940:1:1967;
n_years=size(all_years,2);

%% Useful parameters
n=7^4;
p=n*(n-1)/2;
m_R1=p;
random_indices_pairs_original_R1 = randperm(p, m_R1).';

n=13^4;
p=n*(n-1)/2;
m_R2=10^6;
random_indices_pairs_original_R2 = randperm(p, m_R2).';

n=7^4;
p=n*(n-1)/2;
m_restricted_all_R1=[10^5 m_R1]; 
random_indices_pairs_all_R1=cell(size(m_restricted_all_R1,2),1);
random_indices_pairs_all_R1{end}=random_indices_pairs_original_R1;
for t=1:size(m_restricted_all_R1,2)-1
    random_indices_pairs_all_R1{t}=randperm(p, m_restricted_all_R1(t)).';
end

n=13^4;
p=n*(n-1)/2;
m_restricted_all_R2=[10^5 m_R2]; 
random_indices_pairs_all_R2=cell(size(m_restricted_all_R2,2),1);
random_indices_pairs_all_R2{end}=random_indices_pairs_original_R2;
for t=1:size(m_restricted_all_R2,2)-1
    random_indices_pairs_all_R2{t}=randperm(p, m_restricted_all_R2(t)).';
end

clear random_indices_pairs_original_R1 random_indices_pairs_original_R2

Ueq=0;
Uineq=0;
tol=10^(-4);

param_MOSEK.MSK_IPAR_LOG = 0; 
param_MOSEK.MSK_IPAR_INTPNT_BASIS = 'MSK_BI_NEVER';
param_MOSEK.MSK_IPAR_NUM_THREADS = 1;

%% Construct grid
u0grid=0;
u1grid=1; %it will be replaced in f_s
u2grid=linspace(G_round_larger{2,x}(1),G_round_larger{2,x}(2),100); 
u3grid=linspace(G_round_larger{3,x}(1),G_round_larger{3,x}(2),100); 
u4grid=linspace(G_round_larger{4,x}(1),G_round_larger{4,x}(2),100);
% try different sizes of grid hypercube as discussed in Appendix C of the paper (need powerful cluster to run the code for larger sizes of grid hypercube)
[ca, cb, cc,cd, ce] = ndgrid(u0grid, u1grid, u2grid, u3grid, u4grid);
u0grid=ca(:);
u1grid=cb(:);
u2grid=cc(:);
u3grid=cd(:);
u4grid=ce(:);


%% Workers
workers=600; %number of parallel workers
jobs=round(size(u2grid,1)/workers); %number of jobs per parallel worker

